Cite this Article
Zhao Zhu, Hui Mei, Xia Zheng-Zheng, Zhao Yue-Jin. Centroids analysis for circle of confusion in reverse Hartmann test. Chinese Physics B, 2017, 26(5): 050204  
Permissions
Centroids analysis for circle of confusion in reverse Hartmann test
Zhao Zhu1, 2, Hui Mei1, 2, †, Xia Zheng-Zheng1, 2, Zhao Yue-Jin1, 2
Beijing Key Laboratory for Precision Optoelectronic Measurement Instrument and Technology, Beijing 100081, China
School of Optoelectronics, Beijing Institute of Technology, Beijing 100081, China

 

† Corresponding author. E-mail: huim@bit.edu.cn

Abstract

The point spread function (PSF) is investigated in order to study the centroids algorithm in a reverse Hartmann test (RHT) system. Instead of the diffractive Airy disk in previous researches, the intensity of PSF behaves as a circle of confusion (CoC) and is evaluated in terms of the Lommel function in this paper. The fitting of a single spot with the Gaussian profile to identify its centroid forms the basis of the proposed centroid algorithm. In the implementation process, gray compensation is performed to obtain an intensity distribution in the form of a two-dimensional (2D) Gauss function while the center of the peak is derived as a centroid value. The segmental fringe is also fitted row by row with the one-dimensional (1D) Gauss function and reconstituted by averaged parameter values. The condition used for the proposed method is determined by the strength of linear dependence evaluated by Pearson’s correlation coefficient between profiles of Airy disk and CoC. The accuracies of CoC fitting and centroid computation are theoretically and experimentally demonstrated by simulation and RHTs. The simulation results show that when the correlation coefficient value is more than 0.9999, the proposed centroid algorithm reduces the root-mean-square error (RMSE) by nearly one order of magnitude, thus achieving an accuracy of pixel or better performance in experiment. In addition, the 2D and 1D Gaussian fittings for the segmental fringe achieve almost the same centroid results, which further confirm the feasibility and advantage of the theory and method.

PACS: 02.60.Gf;06.20.-f;07.05.-t;42.30.Va
Keyword:circle of confusion;centroid algorithm;Gaussian approximation;reverse Hartmann test
1. Introduction

Interferometry[1–3] and deflectometry[4,5] have been widely employed in the measurement of optical surfaces. In particular, slope measuring deflectometry (SMD)[6] systems have been developed rapidly. Successful implements include the reverse Hartmann test (RHT) presented by Su et al.[7] and phase measuring deflectometry presented by Knauer et al.[4]

The RHT provides a contact-free, high dynamic range, full field metrology method with easy setup and alignment. It utilizes an incoherent source and imaging detector with external pin-aperture to measure the shape of tested surfaces as shown in Fig. 1(a). In Fig. 1(b), the coaxial design makes the test less sensitive to geometric uncertainty at component position.[8] In RHT, the surface slopes can be calculated by knowing the lit pixel locations ( , on the screen, the pattern centroids ( , on the detector, and the mirror bright region ( , . Since the accuracy of this system geometry is significantly determined by computation of the centroids of ( , , a high-precision centroid algorithm with better reliability becomes very important.

Fig. 1. (color online) (a) Off-axis and (b) on-axis geometry setup of RHTs.

To explore the centroid computing method for the RHT, one needs to start with analyzing the intensity distribution of the point spread function (PSF). The PSF describes the contribution of irradiance on a single image point. In the beam path diagram described in Fig. 1, we use an incoherent imaging system in which beam energy from one object point is imaged as a PSF. In RHT, limited by the diffraction of light, an imaging lens with circular pin-aperture makes the focused spot well described by the Airy disk. In fact, the detector is focused on the reflected surface rather than the generator, which means that the intensity distribution of the PSF onto the image plane behaves as a circle rather than a point. This circle is called the CoC, and it is a measure of how much the image point is defocused. So, we consider the PSF of the imaging detector as a CoC, rather than an Airy disk. In this case, the Huygens–Fresnel integral used to determine the light intensity cannot be obtained in a closed form,[9] and the traditional method of centroid computation[1,7] should be improved.

In this paper, the CoC intensity distribution is determined by the characteristics of imaging lens and pin-aperture model. It is simulated on computer based on Lommel functions and Bessel functions. Image formation experiments are also conducted to prove the correctness and validity of the model. The strength of linear dependence between the Airy disk profile and CoC profile is evaluated by Pearson’s correlation coefficient. Finally, the centroid algorithm of CoC is investigated in the RHT system.

The rest of this paper is organized as follows. In Section 2, the intensity of CoC is analyzed based on the RHT pin-aperture model. In Section 3, a simulation is conducted to verify the performance of the algorithm. In Section 4, the centroid computation experiments for the spot and fringe pattern based on a coaxial RHT system are conducted. In Section 5, some conclusions and a discussion are presented.

2. CoC and Gaussian approximation
2.1. Traditional centroids calculation in RHT

The centroid is generally defined as the irradiance weighted average position of detector pixels. The centroid from pixel-to-pixel light curve data is calculated by[10]

where ( , is the centroid coordinate, ( , and are coordinates and gray of the ith×jth pixel on detector.

Actually, the patterns on the detector are small in size. Figure 2 contains small portions of recorded images, with spots ranging from 7 to 9 pixels in diameter, or with fringes ranging from 14 to 16 pixels in width. The shot noise and dark current noise in the imaging sensor, as well as the stray light from the outside environment will significantly affect the centroid accuracy computed by Eq. (1), and even produce centroid errors.

Fig. 2. Sample patterns recorded on imaging detector in RHT. (a) Enlarged drawing of a spot; (b), (c) Enlarged drawing of partial fringe patterns in x (horizontal) and y (vertical) directions.

To counter the problems above, the Gaussian fitting of the single spot pattern takes into account the sub-pixel localization, and the Gaussian fitting is used to precisely extract its centroid. We know that the Airy disk can be approximated by the Gauss function, however, not all the patterns of RHT can have the same approximation due to diffraction and out-of-focus reasons. Thus we need to establish and evaluate the physical optics model of PSF before the centroid computation in RHT.

2.2. CoC analysis in RHT system

Figure 3 gives an equivalent light path of the coaxial RHT system. An optical generator is used to emit the spots or scan line pattern. This pattern illuminates the reflected surface and the detector captures the tested surface image formed by reflected light rays. The pin-aperture has selective action on passing rays and plays a decisive role in measurement resolution.

Fig. 3. (color online) Light path of coaxial RHT system. d is the pin-aperture diameter, f is the focus in object space, is the focus in image space, is the diameter of CoC.

The light path described in Fig. 3 uses an incoherent imaging system in which energy from one object point is imaged as a PSF. The detector is focused on the reflected surface rather than the generator. It means that the PSF onto the image plane acts as a CoC. Its intensity distribution is determined by the characteristics of imaging lens and pin-aperture model. In geometrical optics it is assumed that light rays travel along straight-line paths. However, diffraction occurs when the light encounters the pin-aperture. Thus, we aim at determining the intensity distribution within the CoC for defocused points.

The pin-aperture system consists of a circular opening in a planar opaque screen. The image plane is placed at a distance or behind the pin-aperture, where is the focal length and is the defocal distance. Furthermore, we assume that the transfer function of the lens is simply a phase transformation and does not affect the transmittance. The complex function of position is given by

where and are the amplitude and phase of the wave at position , respectively. Based on the Huygens–Fresnel principle, the propagation of the wave over the distance , whose extent is limited by the pin-aperture, can be described with the help of a linear space-invariant system. Thus the field amplitude on the detector image plane ( , can be expressed by a two-dimensional (2D) convolution of the input wave field distribution and the impulse response of the system evaluated over the area of the pin-aperture as follows:
where the wave number , is a unit vector in the direction of the optical axis, and is a vector from the center of the aperture to the observation point on the image plane. The is the radial distance in the pin-aperture plane. When , C = 1, otherwise, C = 0. If the distance is large compared with the diameter of the aperture and the region of observation is limited to an area around the optical axis (the Fresnel approximation), then the obliquity factor and reduces to the Fourier transformation of the portion of field subtended by the lens aperture. Assuming , and the light intensity (infinite time-average of the amplitude) for monochromatic light is given in this case by , then
which is the well-known Airy pattern, where J1 is the Bessel function of the first kind and the first order.

In classical physics, the diffraction phenomenon is described as the interference of waves according to the Huygens–Fresnel principle. However, the Huygens–Fresnel integral used to determine the light intensity in the CoC cannot be obtained in a closed form and must be evaluated in terms of Lommel functions and Bessel functions. The light intensity in the CoC is given as follows:

where I0 is the maximum intensity of the pattern, while U1 and U2 are Lommel functions. The Lommel functions can be evaluated by a series of approximations given by
where is the Bessel function of the order, and is expressed as follows:
The variables u and v specify the position of a point within the CoC, and they can be expressed as
where and is the wavenumber. Especially, at u = 0, the distribution of light intensity of a point focused on the image plane is given by the following equation which is the derived Airy pattern:

Simulations of the intensity distribution are plotted in Fig. 4, which shows that the distribution changes significantly with u. When , the intensity distributions become complicated and will cause potential problems in centroid data calculation, so this situation should be avoided in an actual test.

Fig. 4. (color online) Simulations of intensity distribution . Each simulation contains three-dimensional (3D) energy distribution surfaces and CoC patterns for parameter values (a) u = 0 (Airy disk), (b) , (c) , (d) , (e) , and (f) .

An important observation is that when u ranges from 0 to , the CoC keeps a similar profile of Airy disk (Fig. 5). In this paper, the strength of linear dependence between Airy disk profile and CoC profile is evaluated by Pearson’s correlation coefficient. The classical estimator of the correlation coefficient is given by the sample correlation coefficient, commonly denoted as r.

where ( , , ( , , …, ( , are the observed samples of Airy disk profile and CoC profile; and are the average values of sample respectively. The r values range from −1.0 to 1.0. The coefficient of 1.0 indicates a perfect positive correlation.[11] The correlation coefficients between and are listed in Table 1. The r values greater than 0.9999 indicate that there are strong positive relationships between the and when , which accords with the conclusion of Fig. 5.

Fig. 5. (color online) Profile curves of Airy disk and CoCs for parameter values , , , and .
Table 1.

Parameters of the scanning patterns collecting system.

.
2.3. CoC centroiding algorithm with Gaussian approximation profile

As described in subSection 2.1, the light source is blurred in a characteristic manner by the RHT system with a pin-aperture, termed the CoC. When the value of u is kept in a range from 0 to , the point of light appears as a bright maximum cone of light surrounded by concentric diffraction rings in the image plane. The shape of this maximum cone which is similar to the shape of Airy disk can be approximated by a Gaussian fit. This idea from Airy disk approximation using a 2D Gaussian function describes the intensity distribution produced by a point source.[12] The approximation of the PSF by a Gaussian function is relatively simple, computationally quick, and more robust under the condition of a lower signal-to-noise ratio than that from the mass center expressed by Eq. (1).[13,14] It is important to note that computationally a Gaussian fitting can be made by either least squares fitting or maximum likelihood estimation. Maximum likelihood estimation is a more accurate method, but more complex than least-squares fitting.[14,15]

The fitting of a single CoC spot with a Gaussian function to identify its centroid forms the basis of the proposed centroiding algorithm. The Gauss function used to formulate the gray distribution of the spot pattern is given by

where a is the height of the peak, ( , is the position of the center of the peak, which also serves as the calculated centroid, and c1 and c2 are the Gaussian RMS widths.

The algorithm of Eq. (11) is also applied to the fringe pattern if we characterize the fringe as a superposition of multiple continuous spot patterns. In fact, the curve of the 2D Gauss function has a characteristic symmetric “bell” shape, thus the fitting with the 2D Gaussian function about the fringe pattern will present an “oblate bell” shape. In order to further test the accuracy of the fringe centroid, 1D Gaussian fitting is conducted row by row on a segmental fringe, and then the average values of parameters are used to build the Gauss function. In this way, the plot of fitted fringe is a true ‘line’. The 1D Gauss functions used to formulate the gray distributions of fringe patterns are given by

where mean(*) represents averaging operation, and are Gauss functions for fitting the fringe patterns in the vertical and horizontal directions, and are the heights of the peaks, while and are Gaussian RMS widths.

3. Centroid algorithm simulation

A simulation experiment is carried out to test the precision of the proposed centroid algorithm. A simulated spot image is generated by the computer; each spot is defended as a CoC as shown in Fig. 6(a). It contains 10×10 spots, whose centroids are known and served as the truth-values, denoted by (x0, . Different centroid algorithms are used to calculate the centroids, while the noises are artificially added into this simulated image as shown in Figs. 6(b) and 6(c). Table 2 shows the RMSE results of the simulation experiments, in which ( , are the centroids computed by Eq. (1), and ( , are the centroids computed by Eq. (11). It shows that the centroid algorithm with Gaussian approximation can increase the precision by almost one order of magnitude.

Fig. 6. Simulations of centroid computation, showing (a) the standard CoCs image, whose centroids are determined by computer program, (b) the centroid computation result from the method with Eq. (1), and (c) the centroid computation result from the method with Eq. (11). The small blue circles represent the known centroid coordinates. The red crosshairs refer to the centroid coordinates computed by (b) Eq. (1) and (c) Eq. (11).
Table 2.

Parameter of the scanning patterns collecting system.

.
4. Experiments and results
4.1. Layout of experiments

We designed a coaxial RHT to capture target images. The experimental setup consisted of an image generator (light source), PBS, reflector surfaces (planar, spherical, and off-axis parabolic (OAP) mirror), pin-aperture, imaging detector, and a computer as shown in Fig. 7. The parameters are shown in Table 3. The LCD screen was employed as an image generator, and it generated a spot matrix or multiple fringe pattern, whose width and space could be controlled by computer. Then the light from the LCD is reflected by the PBS and reflector surface. After passing the PBS and pin-aperture, rays were imaged on the detector’s receiving surface, which were sent to an image processing system and displayed on the monitor. The recorded patterns were used to conduct the centroid calculation. The coaxial design makes the RHT less sensitive to geometric uncertainty in component positions.[8]

Fig. 7. (color online) Coaxial RHT setup, where the detector captures the image of the reflector and records the spots or fringes pattern.
Table 3.

Parameter of the RHT system.

.
4.2. Experiments for spot pattern
4.2.1. Spots intensity analysis and the Gaussian fitting

The spots patterns with different shapes are collected by the setup shown in Fig. 7. In Fig. 8, the patterns are successively formed by a flat crystal, spherical mirror, and OAP mirror. The spot shape in Fig. 8(a) is circular. However, the ones in Figs. 8(b) and 8(c) are changed into non-circle shapes because of the aberration effect caused by the spherical or aspherical mirror.

Fig. 8. Spot array images for centroid computing in the coaxial RHT experiments. They are formed by (a) flat crystal, (b) spherical mirror, and (c) OAP mirror.

Through the binarization and the statistics for connected regions, the invalid data with small area are removed. Once the valid images are acquired, we perform a 2D Gaussian fitting by least squares using a custom Matlab (The Mathworks, Natick, MA) script. Each circular spot intensity is approximated by a Gaussian fit and compared with the simulation by Lommel function (Eq. (5)). Figures 9 and 10 show the examples of fitting results by profile curves and intensity distribution, which can be quantitatively evaluated by the calculated values in Table 4. It indicates that the original spot pattern contains valid centroids, the root-mean-square errors (RMSEs) of fitting are in a range of about 2.5–2.6 in grey-scale value.

Fig. 9. (color online) Profile curves of the fitting results, showing (a) profiles of a circular spot and profiles of a non-circular spot in (b) x and (c) y cross sections.
Fig. 10. (color online) 2D Gaussian fitting results of spots formed by (a) flat crystal, (b) spherical mirror, and (c) OAP mirror. I0 denotes the initial intensity, and , the intensity fitted by 2D Gaussian function.
Table 4.

Numerical analysis for spot Gauss fit. The calculations are completed by Matlab procedure.

.
4.2.2. Spot centroid computation

The spot centroids are computed by the methods described by Eq. (1) and Eq. (11). Table 5 shows the results in detail, where ( , is the centroid calculated by Eq. (1), and ( , is the difference between ( , and ( , . It can be found that the centroid difference ( , is less than 0.1 pixel. The data from each spot are an average of 25 single images with a total exposure time of 68 seconds per step.

Table 5.

An example of centroid computation results (spots patterns, units: pixels).

.

In fact, the spot reflected by OAP demonstrates oblate shape ( , ). According to Table 5, the Gaussian fitting gains stabilized compensation and centroid calculation results, which proves that the centroid algorithm with 2D Gaussian approximation is applicable to a spot with an oblate shape.

4.3. Experiments for fringe pattern
4.3.1. Fringes intensity analysis and the Gaussian fitting

The fringe patterns are captured and analyzed in the same way as spot patterns. As shown in Fig. 11, the fringes reflected by flat crystal and spherical surface retain their linear characters. However, the OAP surface makes the shape of fringes more complex than the other ones, which are presented as ‘curves’.

Fig. 11. Fringe images for centroid computation in the coaxial RHT experiments. They are formed by (a) flat crystal, (b) spherical mirror, and (c) OAP mirror.

Invalid data are located at both ends of the fringe. In order to improve the efficiency, instead of computing the centroid row-by-row, we segment the successive fringe into several separate short lines and calculate the centroid with 2D Gaussian fitting for each of them. The width of the fringe yields the length of each segment in our research. Figures 12–14 show the examples of segmental ‘straight-lines’ and the fitting results. The numerical computation results in Table 6 show that the RMSEs of fitting are in a range of about 3–7 in grey-scale value.

Fig. 12. (color online) Gray compensation with 2D Gaussian fitting for the fringes in the (a) horizontal and (b) vertical directions. The reflector surface is flat.
Fig. 13. (color online) Gray compensation with 2D Gaussian fitting for the fringes in the (a) horizontal and (b) vertical directions. The reflector surface is spherical.
Fig. 14. (color online) Gray compensation with 2D Gaussian fitting for the fringes in the (a) horizontal and (b) vertical directions. The OAP is used as the reflector.
Table 6.

Numerical analyses for fringe Gauss fit. The calculations are completed by Matlab procedure.

.
4.3.2. Fringe centroid computation

The fringe centroids are computed by the methods described in Eqs. (1) and (11). Tables 7 and 8 show the results in detail. The centroid difference ( , is less than (0.05, 0.07) pixel.

Table 7.

An example of centroid computation results (horizontal fringes, units: pixels).

.
Table 8.

An example of centroid computation results (vertical fringes, units: pixels).

.

In fact, the curve of the 2D Gauss function has a characteristic symmetric ‘bell’ shape. The fitting result about the fringe pattern presents an ‘oblate bell’ shape. Taking Fig. 12(b) for example, the full width at half maximum (FWTM) in the x direction is pixels, in the y direction is pixels. In order to further test the accuracy of the fringe centroid, 1D Gaussian fitting (Eq. (12)) is conducted row by row, and then the average values of parameters are used to build the Gauss function of the fringe pattern. By this method, the shape of fitted fringe presents a true “line”, rather than an ‘oblate bell’. The shape of ‘oblate bell’ arises from the differences between 2D Gaussian fitting and 1D mean Gaussian fitting as shown in Fig. 15. However, the values of are almost identical, which proves that the accuracy of fringe centroids calculated by 2D Gaussian fitting satisfies the requirement.

Fig. 15. (color online) Comparison between (a) 2D Gaussian fitting and (b) 1D mean Gaussian fitting for Fig. 12(b), and (c) the difference between panels (a) and (b).
4.4. Anti-noise experiments

Anti-noise experiments for centroid algorithms with the methods of Eq. (1) and Eq. (11) are conducted. In this test, an initial pattern is captured by a coaxial RHT system with features equivalent to a standard CoCs image, Fig. 6(a) in simulation, and the results of centroid computations with different methods serve as the standard values. Under the same system, by reducing the required exposure times or the gains of the detector, images with different SNRs are obtained, and the results of the centroid computation are used as tested values. The standard and tested values are compared with each other to analyze the anti-noise capabilities of the two centroid algorithms.

Table 9 shows the results of anti-noise experiments , where ( , and ( , are centroids computed by Eqs. (1) and (11), ( , and ( , are centroids shifting between the initial image and noised images computed by these two different methods, ( , is the difference between ( , and ( , .

Table 9.

An example of result in noise-suppression experiment.

.

It can be found that a greater shift of centroid is obtained when the SNR value the noised image has is larger. The value of ( , shows that the centroid shift computed by Gaussian fitting is less than that computed by Eq. (1). The values of ( , are also enlarged with the decrease of SNR, which implies that the centroid shift computed by Gaussian fitting grows at a slower pace relative to Eq. (1). Anti-noise tests prove that the Gaussian fitting method has the better anti-noise capability in centroid computation, which is consistent with simulation result.

5. Conclusions and perspectives

A centroid algorithm with Gaussian approximation is theoretically and experimentally demonstrated based on the study of CoC in the RHT system. Since the detector is not focused on the light source, the intensity of PSF is evaluated in terms of Lommel functions and Bessel functions. According to the CoC computational simulation, we find that when the correlation coefficient r is more than 0.9999, there is a strong positive relationship between the and , which indicates that the centroid can be calculated by using the Gauss profile approximation method. A simulation experiment is conducted to verify the accuracy of the proposed centroid algorithm. The simulation result shows that the RMSEs are reduced by almost one order of magnitude compared with those from the traditional method. The further actual tests on spot and fringe patterns are carried out by a coaxial RHT system. The anti-noise tests further prove that the Gaussian fitting method has the better anti-noise capability in centroid computation, which is consistent with the simulation result. The results of the centroid achieve an accuracy of 0.01–0.1 pixel or better performance. The contrast experiments of 2D and 1D Gaussian fittings for the segmental fringe obtain the same centroid results, which prove that the centroid computation by using the 2D Gaussian method not only has the advantage of high efficiency, but also meets the precision requirements.

Reference
[1] Greivenkamp J E Gappinger R O 2004 Appl. Opt. 43 5143
[2] Baer G Schindler J Pruss C Siepmann J Osten W 2014 Int. J. Optomechatronics 8 242
[3] Zhang S Tan Y D Zhang S L 2014 Chin. Phys. 23 114202
[4] Knauer M C Kaminski J Hausler G 2004 Photonics Europe September 10, 2004 Strasbourg, France 366 10.1117/12.545704
[5] Zhou T Chen K Wei H Li Y 2016 Appl. Opt. 55 2760
[6] Su T Maldonado A Su P Burge J H 2015 Appl. Opt. 54 2981
[7] Su P Parks R E Wang L Angel R P Burge J H 2010 Appl. Opt. 49 4404
[8] Huang R Su P Horne T Brusa G Burge J H 2014 Opt. Eng. 53 085106
[9] Potmesil M Chakravarty I 1981 Comput. Graph. 15 297
[10] Choi J N Ryu D Kim S W Kim D W Su P Huang R Kim Y S Yang H S 2015 Adv. Space Res. 56 2483
[11] Shevlyakov G Smirnov P 2016 Aust. J. Stat. 40 147
[12] Putman C A De Grooth B G Van Hulst N F Greve J 1992 J. Appl. Phys. 72 6
[13] Verdaasdonk J S Stephens A D Haase J Bloom K 2014 J. Cell. Physiol. 229 132
[14] Stallinga S Rieger B 2010 Opt. Express 18 24461
[15] Laurence T A Chromy B A 2010 Nat. Methods 7 338